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Plane Poiseuille flow past a nanoscale cylinder that is arbitrarily confined (i.e., symmetrically or 
asymmetrically confined) in a slit channel is studied via hydrodynamic lubrication theory and molec¬ 
ular dynamics simulations, considering cases where the cylinder remains static or undergoes ther¬ 
mal motion. Lubrication theory predictions for the drag force and volumetric flow rate are in close 
agreement with molecular dynamics simulations of flows having molecularly thin lubrication gaps, 
despite the presence of significant structural forces induced by the crystalline structure of the mod¬ 
eled solid. While the maximum drag force is observed in symmetric confinement, i.e., when the 
cylinder is equidistant from both channel walls, the drag decays significantly as the cylinder moves 
away from the channel centerline and approaches a wall. Hence, significant reductions in the mean 
drag force on the cylinder and hydraulic resistance of the channel can be observed when thermal 
motion induces random off-center displacements. Analytical expressions and numerical results in 
this work provide useful insights into the hydrodynamics of colloidal solids and macromolecules in 
confinement. 


I. INTRODUCTION 

Understanding hydrodynamic phenomena at the nanoscale has become i ncre asingly important with the 
advent of numerous new technologies enabled by nanofabrication met ho d^^^. Predicting hydrodynamic 
forces and volumetric rates in nanoscale flows is particularly relevant to the design of nano-electromechanical 
systems (NEMS), nanowire nanosensors, and nanofluidic devices for applications that range from bioengi¬ 
neering to materials science and renewable energjf^HH, At length scales in the order of one nanometer, which 
corresponds roughly to the size of three water molecules, fundamental assumptions upon which classical 
hydrodynamic equations are derived need to be thoroughly evaluated. For example, experimental studies 
indicate that solid-like transitions and viscoelastic behavior can arise for simple liquids confined in nanoscale 
lubrication gap^SEI] Transitions to viscoelastic behavior of simple fluids are found to occur when hydrody¬ 
namic time scales are comparable to the relaxation time of the fluid, which in micro/nanoscale confinement 
can exceed by several orders of magnitude its bulk valu^^^Ell Another fundamental issue for continuum- 
based descriptions of micro/nanoscale flows is the difficulty to determine proper boundary conditions (e.g., no 
slip, Navier/Maxwell slip) that are determined by physicochemical properties of the fluid and the nanoscale 
structure of the confining solid surfaceii^lti^. In this context, molecular dynamics (MD) simulations have 
become a valuable tool to study hydrodynamic phenomena in nanoscale confinement and to help in devel¬ 
oping and validating continuum-based descriptions for nanoscale flows. Previous works have demonstrated 
the use of different MD techniques (e.g., equilibrium/non-equilibrium) to determine hydrodynamic forces, 
boundary conditions, and transport coefficients resulting from complex interfacial phenomena in nanoscale 
confinemenlpHll]^ Along similar lines, the present work resorts to fully atomistic non-equilibrium MD simu¬ 
lations in order to develop (continuum-based) hydrodynamic descriptions that yield quantitative predictions 
for nanoscale lubrication flows past micro/nanoscale bodies that are perfectly static or subject to thermal 
motion. 

Hydrodynamic lubrication theory is suitable to study flows within lubrication films having nonuniform 
thickness and arbitrary shape. Hydrodynamic models for nanoscale films, however, can encounter signifi¬ 
cant limitations due to complex interfacial phenomena that arise near fluid-solid interfaces. For example, 
nanoscale roughness of the confining solid surfaces can produce lubrication forces resulting from both liquid- 
solid and solid-solid friction. Moreover, layering of liquid atoms at the surface of a crystalline solid and 
the dynamic rearrangement of the induced solid-like structures can lead to strong structural forces and 
so-called “stick-slip” motiorl^^HlSI, 

Numerous experimental studies using a surface force apparatus or atomic 
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force microscope have probed and characterized lubrication forces in molecularly-thin films at different shear 
rate^22H3ll Notably, in the case of liquids with simple molecular structure (e.g., molecular chains such as n- 
alkanes) confined by molecularly smooth surfaces, tribological studies report that the shear viscosity within 
the lubrication film does not vary significantly with respect to th e liquid bulk value for shear rates as high as 
10® 1/s and lubrication gaps as small as ten molecular diametergSMM]^ Furthermore, these studies indicate 
that the slip plane lies within one molecular diameter from the solid surface independently of the presence 
of electrostatic double-layers or structural forces. For the particular case of water confined by smooth silica 
surfaces the same conclusions hold for lubrication gaps as small as 2 nm (i.e., about the size of six wa¬ 
ter molecules )21221t391 Hence, experimental evidence indicates that hydrodynamic lubrication theory is able 
to yield reliable predictions for atomically smooth surfaces and molecularly-thin lubrication gaps (i.e., as 
thin as five to ten molecular layers) under a wide range of flow conditions (e.g., Couette-type flows with 
moderate-to-high shear rates). 

In this work we study creeping flow of a simple molecular liquid past a colloidal solid cylinder confined in 
a slit channel and lying at arbitrary distances from the channel centerline. In Sec. II we first we study the 
case where the cylinder is perfectly static, employing hydrodynamic lubrication theory we obtain analytical 
expressions for the drag forces and flow rates in finite channels. The studied problem involves two lubrication 
gaps of variable height that can become molecularly thin as the cylinder approaches contact with a channel 
wall. In Sec. Ill we describe the MD technique employed and simulations performed to provide a microscopic 
description of nanoscale flows, without relying on conventional continuum assumptions. In Sec. IV we assess 
the validity of the hydrodynamic lubrication approach for the case of static cylinders of micro/nanoscale 
dimensions. Drag forces and volumetric flow rates predicted for symmetric and asymmetric confinement, 
in channels with different lengths, are compared against numerical solution of the Navier-Stokes (N-S) 
equations and fully atomistic MD simulations. In Sec. V, employing predictions obtained for static conditions 
we study the case of a colloidal cylinder that performs random displacements induced by thermal motion. 
The approach in this section indicates ways in which hydrodynamic lubrication models can be applied to 
predict mean drag forces and flow rates for confined nanoscale bodies (e.g., nanoparticles, macromoleucles, 
nanowires, nanobeams) that are subject to thermal motion. 


II. POISEUILLE FLOW PAST A STATIC CYLINDER ARBITRARILY CONFINED 


The geometry of the studied flow problem is illustrated in Fig. a circular cylinder of radius R is fully 
confined within a slit channel of height H, width W ^ H, and length L ^ H. Under studied conditions the 
flow is assumed to be steady, two-dimensional, incompressible, isothermal, and Newtonian; the fluid density 
p and shear viscosity p are thus assumed constant. The cylinder center is located at {x = 0,y = yc) and 
thus lies at a vertical distance 6 x H = H/2 — yc from the channel centerline (cf. Fig0. To characterize 
the studied flow we will employ the confinement ratio 



( 1 ) 


and the dimensionless off-center displacement, or asymmetry parameter. 


6 = 


1 _ yc 

2 H' 


( 2 ) 


In addition, the dimensionless channel length I — LjH will be employed to characterize finite length effects 
on the volumetric flow rate and drag force for long but finite channels (Z ^ I). 

Flow in the a;-direction at constant volumetric flow rate (per unit width) Q [m^/s] is sustained by a driving 
force GjW = {pm —Pout)H + pg{HL — here, pin and Pout are the static pressures at the channel inlet 

and outlet, respectively, and pg is a constant body force active on the fluid phase. The channel is assumed 
to be sufficiently long so that a parabolic velocity profile u(±L/2, y) = Umax{^ — with Umax = iQ/2H is 
established at the channel inlet and outlet. The studied conditions correspond to creeping flows with very 
small Reynolds numbers Re = pUmaxH/p I. The drag coefficient is thus defined as 

Xik,S) = ^, (3) 

where D is the drag force per unit width and U = u(—L/2, y^) is the velocity of the “unperturbed” flow 
velocity upstream of the cylinder center. 
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FIG. 1. Plane Poiseuille flow past a confined cylinder with arbitrary off-center displacement (|(5| < (1 — k)/2) and 
high confinement ratio {k = 2R/H > 0.5). Flow at a constant volumetric rate is driven by a pressure differential 
(Ap = Pin — Pout) and/or a constant body force {pg) in the 2 ;-direction. The channel is considered to be sufficiently 
long to develop a parabolic velocity profile U{y) = 17maa:(l — 5^) at the inlet and outlet (x = ±LI2). 


For the case of symmetric confinement (S = 0) and moderate confinement ratio 0.2 k < 0.5 there are 
available expressions for the drag coefficient obtained from approximate solution of the Stokes equations 
with no-slip boundary conditions: 

47r 

^ ylo - (1 + 0.5A:2 -p + Aek^ + Ask») \n k + -h Bik^ + Bek^ + 

with Ao = 0.9156892732, Ai = 0.05464866, Ag = 0.26462967, As = 0.792986, B 2 = 1.26653975, B 4 = 
0.9180433, Bq = 1.877101, and Bs = 4.66549 as derived by FaxerP^. A similar expression for symmetric 
confinement and k < 0.5 has been derived via analytical solution of the Oseen equations by TakaisPf^. For 
the case of asymmetric confinement (S ^ 0), which has received considerably less attention, a perturba¬ 
tive solution of the biharmonic equation performed by Harper!^ has provided an approximate analytical 
expression for X(k,S) that is only valid for fc <C 1 (i.e., for R H). This work is primarily concerned with 
flow configurations with arbitrary off-center displacements of the cylinder, 0 < |(5| < (1 — fc)/2, and high 
confinement ratio, fc —> 1, where a lubrication flow approximation is valid. 


A. Hydrodynamic Lubrication Theory 

For the prediction of drag forces and flow rates via hydrodynamic lubrication theory we will assume 
Newtonian flow regimes and no slip boundary conditions. Results from numerical solution of the full Navier- 
Stokes (NS) equations and MD simulations for different configurations (fc > 0.5, S > 0, I > 5) will be 
compared against the derived analytical expressions. Results from fully atomistic MD simulations will 
assess the validity of the adopted assumptions in the case of nanoscale flows of simple molecular liquids 
confined by wettable surfaces that are atomically smooth. The local height along the channel is h±{x) = H 
for |a:| > R and 


h±{x) = ± ( 5 ^ — \/i?2 _ a;2 fQj. |3»| < 

where the (-I-) and (—) signs correspond to the lubrication gaps above and below the cylinder, respectively. 
While for 5 = 0 both lubrication gaps are equal, either gap fully closes for |5| = (1 — k)/2. In clearing the 
cylinder the flow rate Q splits into Q- = a{k,6)Q, flowing below the cylinder, and Q+ = (1 — a)Q flowing 
above the cylinder. The split factor a can be determined by equating the pressure drops Ap_ = Ap_|_ = 
p{R) —p{—R) across the bottom and top lubrication gaps, which yields the following equation: 

dx dx 

a{k,S)J ^^ = [1-aik,6)] y 


( 6 ) 
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While the split factor in Eq. Intakes the expected value a(k, 0) = 1/2 in symmetric confinement, a{k, i5) —>■ 1 
for (5 —(A: — l)/2 (i.e., when closing the top gap) and a{k, i5) —?► 0 for 5 —(1 — k)/2 (i.e., when closing the 
bottom gap). 

After establishing the flow rates Q± via Eq. it is straightforward to determine the drag coefficient 
X{k,d) using conventional lubrication analysis. The full derivation of the drag coefficient is presented in the 
Appendix, while the main analytical results are summarized in this section. The drag coefficient predicted 
via lubrication theory is 


= - -^{a[/p(fc,(5) - fsik,S)] - (1 - a)/s(fc,-<5)}. 


(1-^2) 

with the flow split factor given by the explicit expression 


a{k, 6) = 


1 


The shape functions in Eqs. are 

3fc2(l/2-5) 


fp{k,S) = 


4 65/2 


TT 

— + atan 


l + fpik,S)/fpik,-S)- 

( k 


V 261/2 


1 k 


and 




TT 

— -\- &tcLIl 


3k^ 1 


862(1/2-6) (1/2-6) 6’ 


k 

261/2 } 


(7) 

( 8 ) 


(9) 


( 10 ) 


Here, the confinement parameter 6 = (y/ — = (1/2 — 6)2 — fc2/4 is introduced for a more compact 

definition of the shape functions fp and fg accounting for pressure and shear drag contributions, respectively. 

In the limit A; —>■ 1 the dominant contribution in Eq. is due to pressure forces that are proportional 
to the lubrication parameter e~ 2 , here e = {H — 2R)/2R = {1 — k)/k is the nondimensional effective gap 
height. Hence, Eq. predicts two limit cases: 

A(fc^ 1,0) = (11) 


for cases of symmetric confinement 6 = 0; and 

X{k l,6max) = 37re“5 


( 12 ) 


for the maximum cylinder displacement \Smax\ = {k — l )/2 where one of the lubrication gap is fully closed. 
Eq. [^recovers the asymptotic behavior A oc in lubrication flows as e —>■ [P2H1I1_ Notably, X{k — 
l,6maa;) = X{k —>■ l,0)/(2-\/2) and for high confinement ratios there is significant reduction in the drag 
coefficient as the cylinder approaches contact with the top or bottom channel wall. 

Depending on the particular application, either the flow rate Q or induced the driving force G (i.e., 
pressure differential and body forces) is prescribed. When the flow rate is prescribed knowing the drag 
coefficient suffices to predict the drag force D = ^UX where U = {3Q/2H){1 — 6^). When the driving force 
is prescribed, however, it is necessary to predict the volumetric flow rate (per unit width) Q{k, 6, 1) in order 
to predict the drag force and the hydraulic resistance for different confinement configurations and channel 
aspect ratios. The lubrication flow approximation yields 

S' i\ {Pin ~ Pout)/L + pg ^ ^ ^ ^ 

Q{k, 6,1 ) = - -^ 

determined by the flow correction factor 


l + a(k,5)fp{k,5)-k 


(14) 


where I = LjH is the dimensionless channel length, and a and fp are given by Eqs. According to 

Eqs. 13 ^ the volumetric rate Qoo for Poiseuille flow is recovered for I —)■ oo and the now rate vanishes 
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Q —>■ 0 for fc —1. It is worth noticing that for long but finite channel lengths (1^1) the flow rate increases 
as the cylinder is displaced from the channel centerline (d > 0) according to the ratio 


Q{k,S,l) _ I + fp{k,0)/2 - k 
Q{k,0,l) I + a(k,6jfp(k,6) — k 


(15) 


A few comments are in order about the expressions derived for the drag coefficient and drag force. For the 
particular case of symmetrically confined cylinders, predictions from Eqs.[7- 10 for the drag coefficient A(fc, 0) 
are in close quantitative agreement with asymptotic formulas for /c —>■ 1 proposed in previous worlP^. For 
asymmetrically confined cylinders, Eq.j^predicts a significant decrease in the drag coefficient X{k, S) as |d| —?> 
(1 — k)/2. Moreover, the derived formulas predict a maximum drag force D for symmetric confinement and 
significant drag reduction in asymmetric confinement with reduction ratios that depend on the dimensionless 
channel length I = L/H. To the best of our knowledge an analytical expressions analogous to Eq. and 
Eq. 13 valid for cylinders in Poiseuille-type flows for high confinement ratio {k > 0.5) and arbitrary off-center 


displacement (0 < |(5| < {k — l)/2), are not available in the previous literature. 


III. MOLECULAR DYNAMICS SIMULATION 


Following standard techniques for non-equilibrium MD simulation^^^I^, 
atoms of species s and s' is governed by a generalized Lennard-Jones (LJ) 


the interaction between any two 
potential 





(16) 


Here, Vij = |ri — rjl is the separation between any two atoms {i,j = 1, N), a is the repulsive core diameter, 
and e is the depth of energy potential minimum, which lies at = (2/Ass'Y^^cr. The simulated system 
is composed of three atomic species that correspond to the fluid (s = 1), cylindrical particle (s = 2), 
and channel walls (s = 3) (cf. Fig. [^. The symmetric attraction coefficients Agg' = control the 
degree of wettability of the modeled solid surfaces and the shear-dependent hydrodynamic slip length; 
the parametrization employed (A12 = A13 = 0.8) produces highly wettable solids exhibiting very small 
hydrodynamic slip on flat or curved surfaces over a wide range of flow condition^^SHSU pgr the simulations 
in this work, fluid atoms conform dimer molecules bound by Finitely Extensible Nonlinear Elastic (EENE) 
potentials 


UpENsirij) 


1 

2 


kpr 


2 

max 


log 



(17) 


where kp is the stiffness of the modeled molecular bond and r^ax adjust its maximum extension. The 
use of FENE potentials in addition to LJ interactions allows further control of rheological properties and 
the volatility of the modeled fluid. A Nose-Hoover thermostat maintains constant fluid and solid atoms at 
temperature T = e/kp {kp is the Boltzmann constant). At initialization the atoms of all species are arranged 
in face-centered cubic (fee) lattice with constant spacing Ax = (see Fig. [^), where n = O.S/cr^ was 

the number density employed in all MD simulations in this work. The mean mass density of the fluid 
p = 0.8m/a^ is constant (here m is the atomic mass), for the modeled conditions the shear viscosity is 
H = 2.9y/W€l<7^. 

Non-equilibrium MD simulations are performed to study drag forces and volumetric rates for Poiseuille- 
type flow past a nanoscale cylindrical particle confined at arbitrary off-center displacements S. The net 
force on the cylindrical particle (s = 2) is F = — ^ dujjj /dx obtained as the sum of all atomic interactions 
with the fluid (s = 1) only; i.e., direct atomic interactions between the cylindrical particle and the channel 
walls are neglected in our MD simulations. Different nanochannels with heights H = 26-90Aa; and lengths 
L = 300-1000Aa: are employed in MD simulations in order to characterize the drag and flow rates in a range 
of confinement ratios (fc = 0.6-0.84) and channel aspect ratios {I ~5-38); a constant width H = lOAx is 
employed in all cases. We consider the solid-liquid interface to be located at the zero isopotential contour 
Ulj — 0 for the solid species (i.e., particle and walls); the confinement ratio k and dimensionless channel 
length I in MD simulations were calculated following this criterion. A constant body force f = mgi is applied 
to each fluid atom in order to drive the flow in the a:-direction and periodic boundary conditions are applied 
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FIG. 2. Geometric setup and flow features in MD simulations, (a) Side and perspective views of the simulation 
domain and range of dimensions employed; the side view shows the initial fee atomic lattice with constant spacing 
Ax = (b) Dimensionless mass density p/p where p = O.Sm/cr® is the mean bulk density; solid-like structure 

and layering of fluid atoms is observed near the solid surfaces, (c) Dimensionless momentum density magnitude 
p|u|/p|uo| where uo is the maximum velocity in the symmetric confinement case (5 = 0). As the bottom gap closes 
5 —^ (1 — fc)/2 = 0.15 the flow through the upper gap becomes twice the value observed in symmetric confinement. 
Reported quantities in (b-c) are obtained via time average and spatial average in the z-direction, for k — 0.7 and 
Z = 5 in cases of symmetric confinement 5 = 0 (left panels) and asymmetric confinement at 5 = 0.11 (right panels). 


in the x and z directions (cf. Fig. [^. In all cases the applied driving force produces flows with low Reynolds 
numbers Re = pUmaxH/ki- ^ 0-3, where Umax = pgH^/8p. Distinctive features of the mean mass and 
momentum density fields simulated via MD are reported in Figs. ^b-c). Near the liquid-solid interface we 
observe layering of fluid atoms near the solid surfaces and a small but finite amount of hydrodynamic slip 
that varies locally. As the cylinder approaches contact with a channel wall (cf. Fig. [^) the flow through 
the narrower lubrication gap rapidly decreases as quantitatively predicted by Eq. 


IV. STATIC CYLINDERS 


Predictions from lubrication theory are compared against numerical sim ulat ions via finite element solution 
of the steady-state N-S equation^^ and MD techniques described in Sec. Ill In our numerical simulations 
periodic boundary conditions are applied at the channel inlet and outlet, a constant body force in the x- 
direction results in a total force magnitude G = pg{HL—'KR?)W driving the flow. While physical conditions 
modeled in (continuum-based) N-S simulations correspond to macroscopic channels where H/a ^ 1, the 
conditions modeled in MD simulations correspond to channels with nanoscale dimensions Hja = 20-90 
(i.e., H ~ 5-30 nm). The reported drag forces D are obtained by subtracting the cylinder weight and 
buoyancy force from the total force computed in numerical simulations and the unperturbed velocity U = 
{?>Q/2H){\ — S^) is determined from the numerically computed flow rate (per unit width) Q. In the case of 
MD simulations, reported quantities correspond to averages over sufficiently long times > 0.2HL/Q for 
which convergence of the reported mean values (within a 10% deviation) is observed after reaching a steady 
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FIG. 3. Drag forces and drag coefficients; theoretical predictions (solid/dashed lines), N-S simulation (open markers), 
and MD simulations (filled markers), (a) Normalized drag force D{k,5,l) — D/jj,{3Qoa/‘2H) = (j)\ for symmetrically 
confined cylinders {5 = 0) versus confinement ratio k = 2R/H for I = L/H = 5-32. The flow correction factor 
cf){k, S, 1) is given by Eq. |14[ Plotted for comparison (in both panels) are predictions for infinitely long channels 
(/ —^ oo and 0=1) based on Faxen’s formula (Eq. and lubrication theory (Eq. [see legend]. For a finite 
driving force G in the limit fc —>■ 1 where Q —>■ 0, the drag force becomes equal to the driving force D = G and 
thus D ^ 81 (dashed lines), (b) Drag coefficient X{k,S) = D/fiU where U = {3Q/2H){1 — 5^) as a function of 
£ = (H — 2R)/2R = (1 — k)/k, for 5 = 0 (i.e., symmetric confinement) and jdj = (1 — k)/2 (i.e., limit case of 
asymmetric confinement where the cylinder contacts either channel wall). For e —>■ 0 the drag coefficient exhibits the 
asymptotic behavior predicted in Eqs. |11[|12| with significant drag reduction in asymmetric confinement. 



FIG. 4. Drag reduction and flow enhancement (i.e., hydraulic resistance reduction) in asymmetric confinement; 
theoretical predictions (solid lines), N-S simulation (open markers), and MD simulations (filled markers), (a) Drag 
coefficient ratio X{k,5)/X{k,0) as a function of the dimensionless off-center displacement 5 for three different con¬ 
finement ratios k. Significant drag reduction is observed in asymmetric confinement |5| > 0 as the confinement 
ratio increases (fc —> 1); as expected, agreement between simulations and lubrication theory predictions from Eq. 
increases for fe —> 1. (b) Flow rate enhancement Q{k, S,l)/Q{k,0,l) as a function of the dimensionless off-center 
displacement <5 for for three different finite channels [see legend]. The hydraulic resistance of the channels increases 
when the confined cylinder moves away from the center. 


flow rate. 

We first analyze the results for the hydrodynamic drag force If on a cylinder in symmetric confinement 
conditions (J = 0). The predicted drag force is D{k,0,l) = /i(3Qoo/2I7)0(fc, 0, Z)A(fc, 0) where A is given 
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by Eq. 0 and (p is given by Eq. here, Qoo = GH^/12fj,LW is the flow rate expected for an infinitely 
long channel {I —^ oo) for the finite driving force G applied in numerical simulations. As showed in Fig. §1, 
theoretical predictions for the drag force as a function of the confinement ratio k in symmetric confinement 
conditions are in close agreement with both numerical solution of the N-S equations and MD simulations 
for long channels with dimensionless length I = L/H = 5-32. For fc —>■ 1, as both (bottom/top) lubrication 
gaps close and the flow vanishes {Q —>■ 0) the force on the cylinder balances the applied driving force 
D{k 1,6,1) = G (cf. [^. It is worth noticing (cf., Fig. [^) that for moderate confinement ratios {k ~ 0.5- 
0.7) as the dimensionless channel length increases {I > 10) the hydrodynamic drag force on the cylinder 
becomes less than half the force applied to drive the flow {D/G < 0.5). 

The drag coefficient X{k, 6 ) = D/ predicted by lubrication theory (Eq.[^ is compared against numerical 
simulations in Fig. Wp for the case of symmetric confinement, where 5 = 0, and the limit case when the 
cylinder contacts either one of the channel walls, where |5| = (1 — fc)/2. Results from N-S and MD simulation 
confirm the expected asymptotic behavior A oc where e = (H — 2R)/2R, for fc —>■ 1 that is predicted 

by Eqs. 11 ■ 12 In the limit case |5| = (1 — fc)/2 where the cylinder contacts one of the channel walls there 
is a reduction of about 65% with respect to the drag coefficient in symmetric confinement (cf. Fig. &)• 
Numerical solution of the N-S equations and MD simulations confirm a gradual reduction in the drag 
coefficient as the dimensionless off-center displacement S increases (cf. Fig. |^,). The drag coefficient ratio 
A(fc, 5)/A(fc, 0) quantifying the reduction of drag in asymmetric confinement as a function of S is reported 
in Fig. for three different confinement ratios fc = 0.6,0.69,0.84. As expected the agreement between 
lubrication theory and numerical simulations improves for large confinement ratios (fc > 0.7). Notably, drag 
coefficients computed from MD simulations of flows having molecularly thin lubrication gaps (i.e., three to 
ten atomic layers thin) are in good agreement with numerical solutions of the full Navier-Stokes equations 
and analytical predictions from lubrication theory adopting no-slip boundary conditions (cf. Fig.|^). 

The lubrication analysis in Sec. |II A| also predicts a gradual reduction in the hydraulic resistance of the 
channel as the off-center displacement of the confined cylinder increases. This effect is observed in simula¬ 
tions as an enhancement in the flow rate Q for a prescribed driving force G as the dimensionless off-center 
displacement 6 increases. As showed in Fig|^, theoretical predictions from Eq. [^for the flow enhancement 
ratio Q{k,S,l)/Q{k,0,l) are in close agreement with numerical simulations for different confinement ratios 
(fc = 0.6, 0.69, 0.84) and finite channels with different lengths (1 = 5.6, 33, 38). Simulations confirm predic¬ 
tions of Eq. |15| the flow enhancement for asymmetrically confined cylinders increases for large confinement 
ratios (fc —>■ 1) and decreases for long channels {I —>■ oo). 


V. COLLOIDAL CYLINDERS 


The lubrication analysis in Sec. |II A produced analytical predictions for the position-dependent drag 
force assuming flow past a perfectly static cylinder. This section discusses the application of hydrodynamic 
lubrication to confined colloidal cylinders that undergo Brownian motion and can be strongly influenced 
by colloidal interactions (e.g., van der Waals attraction, steric repulsion, oscillatory structural forces). The 
employed approach is generally applicable to colloidal particles, whether these are freely convected or bound 
to an equilibrium position by diverse restoring forces. The formulas presented in this section, via invoking 
analytical predictions from Sec. IIA aim to predict the mean (noise-averaged) drag force experienced by 
nanobeams, nanowires, or colloidal probes that constitute a key component of NEMS and nanowire-based 
sensors and actuators. 

For overdamped Brownian motion (i.e., neglecting inertial and memory effects), the mean drag force {D) on 
a confined colloidal particle can be estimated by ensemble averaging the (position-dependent) drag for static 
conditions over the sequence of random displacements induced by thermal motion. Since the drag predicted 
under static conditions varies only in the vertical direction (y-direction), uncorrelated thermal motion in the 
flow direction (x-direction) is not expected to affect the mean drag force. Vertical random displacements can 
be statistically described by a probability density QoiS, t) = g{ 6 , t|5o, to); here So is the initial displacement at 
time to where Qo{ 6 ,to) = 5{6 — 6 o)- Hence, the mean drag force expected for overdamped Brownian motion 
is 


{D{k,l,t)) 


3/x 


r+Sr, 


— / Q{k, 6 ,l)X{k,S)go{S,t)dS, 


i-s„ 


( 18 ) 


where A(fc, 6 ) is the drag coefficient (Eq. derived for static conditions and 6 max 


(1 — fc )/2 as before. 
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Similarly, it is useful to define a noise-averaged drag coefficient 

(A(fc, t)) = r \ik, 6)go{S, t)dS, (19) 

J-Srr.a^ 


in order to characterize the mean drag force {D(k, I, t)) for the case where the flow rate Q is prescribed; this 
case also corresponds to prescribing the driving force in sufficiently long channels {I —t oo) for which the 
flow correction factor (Eq. 141 becomes unity (j){k, S, I —)■ oo) = 1. 

Via MD simulations we analyze the case of a confined colloidal cylinder immersed in a fluid with constant 
thermal energy ksT. A linear restorative force Fg = Ks\/{x — A/2)2 + (y ~ H/2Y is applied to bring 
the colloidal cylinder to equilibrium at the center of the channel where 5 = 0; the “spring” constant is 
varied in the range Kg = 0-1 kgTj in order to modulate the root-mean-square (rms) amplitude of the 


dimensionless off-center displacement Srmsit) = a/(( 5(t) — 5o)^). The mean drag and drag coefficient defined 
in Eqs. [TsHT^ are expected to depend on the dimensionless rms displacement Srms = drms observed in MD 
simulations. In order to produce analytical predictions we will assume that the colloidal cylinder follows an 
Ornstein-Uhlenbeck process and thus 


QoiS,t) = 


^/ ^TT Srms (^) 


exp < - 


5 — 5o exp(t/r) 


1 2 




( 20 ) 


where <o = 0 t = D/Kg is an unknown effective diffusion coefficient. Since boundary effects are 
neglected, the probability in Eq. [^is only valid for small rms displacements Srms{t) ^ Smax = (1 — k)/2. 
Eor fi nite values Kg there is long-time limit t»T where initial condition is forgotten and the probability 
in Eq. 20 becomes g{S,t) = exp[—(5/-\/25j.ms)^]/'\/^5rms- 




0 0.005 0.01 0.015 0.02 0.025 0.03 0.035 0.04 


o 1.05p- 
^ 1 ’'' ' 
^ 0.95 - 


0.9 
^ 0.85 


-k = 0.84 


(C) 


0 0.005 0.01 0.015 0.02 ^ 0.025 0.03 0.035 0.04 



FIG. 5. Drag reduction induced by thermal motion and structural forces for a colloidal cylinder symmetrically 
confined in a slit channel, (a-c) Mean drag coefficient reduction {X{k, Srms)}/^{k, S^ms) versus rms off-center dis¬ 
placement. Plotted lines indicate predictions from Eq. for k — 0.64,0.78,0.84. Filled markers correspond to 
results from MD simulations: (a) (a) k = 0.64, I = 3.56 {xrms=Q)\ (b) (•) k — 0.74, I = 3.5, (►) k — 0.74, I = 3.5 
(a:rms=0), (■) k = 0.78, I = 5; (c) (♦) k = 0.82, I = 5, {j) k = 0.84, I = 5 {xrms=0), (★) k = 0.86, I = 3. (d) Time 
evolution of rms displacement versus dimensionless time t/To {To = gkH^/2kBT) for MD simulations in panel (a); 
solid line is an exponential ht, open markers correspond to values computed from MD simulation, (e) Vertical off- 
center displacement in lattice units Ax versus dimensionless time t/To for MD simulations in panel (a); six different 
realizations showing metastable states, (f) Stationary probability distribution q{S) computed from absolute value of 
displacement-time trace in panel (e) via ensemble average over six realizations, (g) Free energy U{5) computed from 
probability distribution in panel (f). 


Predictions for the mean drag coefficient (Eq. |19[ ) via adopting the probability density in Eq. 
t»T and Srms ^ Smax are compared in Eig. [^against results from MD simulations for colloidal cy 
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symmetrically confined in nanoscale channels of various heights H = 30-90(7. For the MD simulations 
reported in Fig. the cylinder is allowed to “freely” drift in the vertical direction while the motion is 
prescribed in the x-direction; the expected rms vertical displacement is S^ms = 2Z3t (for t > t) and the 
(top/bottom) lubrication gaps in these MD simulations become as small as two atomic diameters. For the 
MD simulations reported in Figs. [^-c a restorative force with different strengths [Kg = 0.1-1 ksT/a^) is 
applied; in this case the maximum rms displacement is bounded, Srms = ^/ksT/Kg, and lubrication gaps 
in this case are always larger than five atomic diameters. In all cases, MD simulations (cf. Fig. report a 
decay in the mean drag force as the rms displacement 6 rms increases, in close agreement with predictions 
from Eq. |19| However, the rms displacements reported in Figs. ^-c as computed from MD simulations 
can become significantly smaller than the expected values for diffusion in homogeneous fluid media when 
the cylinder reaches within three atomic layers from the channel walls. Moreover, for the case where the 
cylinder is free to drift vertically we observe a slow exponential relaxation determined by a diffusive time 
T ^ Tjj = ^kH^/2kBT (see Fig. W)- In fact, when a lubrication gap becomes thinner than five atomic 
layers the displacement-time trace of the cylinder center-of-mass exhibits long-lived metastable states (cf. 
Fig.§i) that indicate the local energy minima that lie at average separation Ax = {p/m) (cf. Fig. [^). 

Assuming Boltzmann statistics, the stationary probability distribution {t —>■ oo) is q{5) = Z~^ exp[—[/( i5)/A:bT], 
where U{5) is the (space-dependent) free energy and Z is the proper normalization constant. A strongly 
non-Gaussian probability distribution computed from the displacement-time trace reveals an oscillatory free 
energy U{6) = —fcBFlog(p) -|- const, that decays away from the walls. This observation explains the poorer 
agreement observed in Fig. between MD simulations and predictions adopting a Gaussian probability in 
Eq. Invalid for free Brownian motion in the long-time limit t ^ Tp. Given that our MD simulations do 
not include atomic interactions between solid atoms, the oscillatory free energy variations are attributed to 
the structural rearrangement of fluid layers caused by the cylinder motion. Hence, the modeled steric and 
van der Waals interactions between solid and fluid atoms induced significant energy barriers (At/ = bkgT) 
and long-lived metastable states when the cylinder is close to the wall. 


VI. CONCLUSIONS 

A hydrodynamic lubrication approach was presented to predict drag forces and volumetric rates for plane 
Poiseuille flow past a confined static cylinder as a function of the confinement ratio fc, the dimensionless 
off-center displacement of the cylinder d, and the dimensionless channel length 1. Analytical expressions for 
the drag coefficient introduced in this work are valid for moderate to large confinement ratios (fc > 0.5) and 
arbitrary off-center displacements (0 < |d| < (1 — k)/2). In the high co nfine ment limit k —>■ 1 the derived 
expressions recover the asymptotic behavior reported in previous workd^^*^. The set of derived formulas 
applies to cases when either the volumetric flow rate or the driving force is prescribed. In addition, the 
derived expressions valid for finite channels are suitable for predicting drag forces and volumetric rates for 
flow past periodic arrays of cylinders. 

As the cylinder moves away from the channel centerline and one of the lubrication gap closes, either above 
or below the cylinder, the flow through the closing gap vanishes and so does its contribution to the drag 
force. The derived expressions quantitatively predict that (i) drag forces and drag coefficients have their 
maximum value in symmetric confinement ((5 = 0), and (ii) there are significant reductions in both the drag 
force and drag coefhcient in asymmetric confinement, as the cylinder approaches either one of the channel 
walls (|(5| —> 1/2 — k/2). Gonversely, the hydraulic resistance for a given confinement ratio k is minimized 
when the cylinder contacts a wall and |(5| = 1/2 — k/ 2 . 

In the case of static cylinders, analytical predictions for the drag force and flow rates are in good agreement 
with numerical solutions of the Navier-Stokes equations and fully-atomistic MD simulations of nanoscale 
channels. Notably, conventional hydrodynamic descriptions adopting no-slip boundary conditions produced 
reliable predictions despite the presence of significant steric effects and structural forces observed in fully 
atomistic simulations with wettable solids. The results in this work indicate that hydrodynamic lubrication 
theory can produce reasonable predictions for molecularly thin lubrication gaps (i.e. down to three atomic 
layers) in the case of Poiseuille-type flows on plane channels with surfaces that are molecularly smooth and 
highly wettable by simple molecular liquids. In fact, MD simulations and continuum models (i.e., lubrication 
theory and numerical solution of the N-S equations) reported comparable values of the drag force when one 
of the lubrication gaps became vanishingly small. The observed agreement, however, can be attributed 
to the fact that the flow through the narrowest lubrication gap decreases, as quantitatively predicted by 
equating the pressure drop through each gap, and the dominant contribution to the drag force comes from 
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the widest lubrication gap, which in all studied cases remained thicker than two atomic layers. The MD 
simulations reported a small hydrodynamic slip that depended on the local surface curvature and shear rate 
magnitude but this effect did not affect significantly the agreement with analytical predictions adopting 
no slip boundary conditions. The presented lubrication analysis can be readily extended to systems with 
partially wettable solids where significant hydrodynamic slip is present, provided that the slip length is a 
known parameter. 

The studied case of nanoscale cylinders undergoing thermal motion revealed a few important effects. For 
symmetrically confined colloidal cylinders the mean (noise-average) drag force, determined via ensemble or 
time average, can be significantly lower than the drag force predicted for a static cylinder. The mechanism 
for the predicted drag reduction is not attributed to hydrodynamic slip but rather to the colloidal cylinder 
randomly moving to off-center positions where the drag predicted in static conditions is significantly lower. 
Similar thermally-induced effects produce a noticeable reduction in the mean hydraulic resistance as the 
rms displacement increases. For creeping flows and after assuming a Gaussian probability density for the 
thermally-induced displacements of the colloidal cylinder, the reduction in the mean drag and hydraulic 
resistance can be quantitatively predicted by averaging the position-dependent drag and flow rate derived 
in static conditions. The observed effects can be enhanced by increasing the rms amplitude of the cylinder 
displacement via different mechanisms, which can include mechanical or acoustic actuation and/or increasing 
the fluid temperature. Under studied conditions, the presence of significant structural forces was found to be 
the major obstacle to safely extending hydrodynamic lubrication theory to nanoscale flows in plane channels. 
Analytical or numerical solution of a Fokker-Planck equation can predict the probability density of random 
thermal displacements but this will require a priori knowledge of local free energy variations for a confined 
colloidal cylinder. Although structural forces did not play a significant role when the cylinder position was 
prescribed, oscillatory structural forces induced strongly non-Gaussian probability densities and long-lived 
metastable positions of the cylinder at integer number of atomic layers from the channel wall. The analysis 
and results presented in this work are relevant to the design of NEMS and nanowire-based sensors and 
actuators, nanofluidic devices for transport and separation of nanoparticles or macromolecules, and can 
potentially guide experimental studies of the nanorheology of confined fluids using colloidal probes. 
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Appendix A: Drag force and flow rate derivation via lubrication theory 


The steady-state, incompressible, and isothermal Navier-Stokes equations in the lubrication limit are 
reduced to: 


d^u dp n n 

oy^ ox ay 


(Al) 


Here, p is the shear viscosity, p is the fluid velocity, and g is acceleration due to a constant body force. 
Solution of Eq. |A1| gives the velocity profile 


(A2) 


and the volumetric flow rate (per unit width) 


Q = 


h^{x) ( dp 


12 p \dx 


pg 


(A3) 


As showed in Fig. the local height is h{x) = H for |a:| > R and h{x) = h±{x) for |a;| < R; here 
h±{x) = ± S) H — ^/R"^ — x"^ where the (-I-) and (—) signs correspond to the top and bottom lubrication 
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gaps, respectively. In clearing the cylinder the flow must split into Q- = a(k,5)Q below the cylinder and 
Q+ = a{k, —S)Q above the cylinder; the split factor a must satisfy mass conservation and thus a(k, —S) = 
1 — a{k,S). The pressure gradient inside the bottom lubrication gap is given by 


dp- 

dx 


pg - l2pa{k,5)Qj^ 


(A4) 


Exact analytical integration of Eq. A4 gives the pressure drop across the bottom gap 


Ap_ = p{R) - p{-R) = pg2R - 


l2pQ 


o:{k,S)fp{k,d) 


(A5) 


where 


fpik,S) = - 


3k^{l/2-S) 


55/2 


TT f k 


3fc3 


+ 


8&2 {1/2-S) {1/2-5) b 


(A6) 


and b = (1/2 — 5)'^ — fc^/4. Similarly, the pr essur e drop across the top gap is Ap+ = pgHk — [1 — 
a{k,d)]fp{k,—6). In the lubrication limit (Eq. |A1[ ) the pressure drop across the bottom and top gaps 
must be equal, Ap_ = Ap+ = Ap, and thus the now split factor is given by 


a{k, 5) = 


1 


l + fp{k,6)/fp{k,-S) 


Using Eq. A7 the velocity profile for |a;| < i? in the top/bottom (+/—) gap can be cast as 

u±{x,y) = a{k,T5)-^{h±y - y"^), 


(A7) 


(AS) 


and thus the shear stress on the (top/bottom) channel wall is 


T±{x,y) = =Fa(fc, T<5) 


&p.Q 


(A9) 


Analytical integration of Eq. A9 wall segments above and below the cylinder {—R < x < R, y = 0 and 
y = H) gives a shear force (per unit width in the a;-direction) 


where 


Fs = [a{k,5)fs{k,5) 


oi{k,-5)fs{k,-5)] 


fs{k,S) 


4b3/2 


TT 

— + atan 



k 


(AlO) 


(All) 


To find the drag force D (per unit width) on the cylinder for a given flow rate Q we apply a control 
volume approach. Static force equilibrium in the x-direction within the channel section containing the 
cylinder(—i? < x < R) gives 

[p{-R) - p{R)]H + pg{2RH - ttR^) + F, + F = 0 (A12) 


where F = —{D + Fi,) is the force exerted on the fluid by the cylinder and Ff, = —pgirR/^ is the buoyancy 
force (per unit width). Introducing Eq. A5 and Eq. AIO into Eq. A12 gives 


D = 


12pQ 

H 


{a{k, 5) [fp{k, (5) - fs{k, <5)] + [I - a{k, 5)] fs{k, -<5)} . 


(A13) 


This analysis further considers the case where the flow rate Q is not prescribed but rather determined 
by pressure differentials and body forces applied. When the obstruction to the flow due to the confined 
cylinder is negligible, which corresponds to the limit of infinitely long channels, the flow rate is Qoo = 
[{Pin — Pout)/L + pg\H^/12p as predicted for plane Poiseuille flow. The additional hydraulic resistance 
caused by the cylinder produces a reduction of the flow expected for plane Poiseuille flow and thus the 
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actual flow rate is Q = (j){k,6,l)Qao: where (/> < 1 is a flow correction factor that needs to be calculated 
for different confinement configurations and channel aspect ratios. Assuming parabolic velocity profiles are 
recovered at the channel ends, static equilibrium for the fluid contained in the entire channel gives 


[p^n-Pout)H + pg{2LH - 7 tR^) + F, - _ 2i?) + A = 0, 


(A14) 


where pin and Pout are the prescribed pressures at the channel inlet and outlet, respectively. Combining 
Eq. A12 with Eq. A14 gives the flow rate 


Q{k, S, 1) = 


ijpin - Pout)H^ + pg2LH‘^ 

Vlp, [I + a{k, S)fp(k, S) — k] 


(A15) 


where I = L/H is the dimensionless channel length or longitudinal aspect ratio. Using Eq. A16 for the flow 
rate we can define the flow correction factor 


ms,i) = 


Q{k,S,l) 

Qoq 


I 


[I + a{k,S)fp{k,S) - k]’ 


(A16) 


as expected for channels that are not fully blocked by the cylinder (fc < f) we have 0 —)■ 1 in the limit 

I —>■ 00 . 
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